\documentclass[12pt]{amsart}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{graphicx}

\usepackage[ruled]{algorithm2e}
\SetKwInput{KwInit}{Initialize}

\usepackage{multicol}

\usepackage{breqn}

%make things fit on less pages:
\usepackage{times}
\usepackage{fullpage}
\usepackage[small,it]{caption}
%\usepackage{savetrees}
%\usepackage[small,compact]{titlesec}
\usepackage{setspace}
%\addtolength{\belowcaptionskip}{-5mm}
\addtolength{\intextsep}{-5mm}


\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}

\newtheorem{thm}{Theorem}[section]
\newtheorem{lem}[thm]{Lemma}

\theoremstyle{remark}
\newtheorem{rem}{Remark}
\newtheorem{defn}{Definition}


\title{Hybrid regularization for MRI reconstruction with field inhomogeneity correction}
\author{Ryan Compton, Louis Bouchard, Stan Osher}
\date{October 96, 2010}


\begin{document}

\begin{abstract}
Rapid acquisition of MR images via reconstruction from undersampled $k$-space data has the potential to greatly decrease MRI scan time on existing medical hardware. To this end iterative image reconstrction based on the technique of compressed sensing has become the method choice for many researchers \cite{Lustig2007}. However, while conventional compressed sensing relies on random measurements from a discrete Fourier transform, actual MR scans often suffer from off-resonance effects and thus generate data by way of a non-Fourier operator, complicating image reconstruction methods and introducing additional computational bottlenecks \cite{Fessler2005}.

In this work we introduce a combined TV-framelet regularization to the undersampled MR reconstruction problem in inhomogeneous fields. We show that the inclusion of a framelet regularization term decreases computation time and improves image quality over the traditional total variation based approach.

\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
Our model of MR image formation is based on solving for the proton density map, $\rho(\mathbf{r})$, in the signal equation

\begin{equation}\label{sigeq}
s(t) = \int_{\mathbb{R}^2} \rho(\mathbf{r})e^{-z(\mathbf{r})t}e^{-2\pi i \mathbf{k}(t) \cdot \mathbf{r}} d\mathbf{r}
\end{equation}
where $s(t)$ is the recorded signal, $z(\mathbf{r})$ contains information about relaxation and off-resonance effects, and the $k$-space coordinates, $\mathbf{k}(t)$, are known only on a nonuniform and possibly undersampled grid \cite{Haacke1999}. An $n$ point discretization of the physical space integral at $m$ different values of $t$ leads to the $m \times n$ linear system

\begin{equation}\label{dissigeq}
\vec{s} = A \vec{\rho}
\end{equation}
with system matrix
\begin{equation}\label{aij}
A_{ij} = e^{-z(\mathbf{r}_j)t_i}e^{-2\pi i \mathbf{k}(t_i) \cdot \mathbf{r}_j}.
\end{equation}
The goal of MRI reconstruction is to efficiently invert the linear system (\ref{dissigeq}) in a way that produces high quality images.

Classically, the background field is assumed perfectly homogeneous and the phase accrual due to off resonant frequency, $e^{-z(\mathbf{r}_j)t_i}$, is then ignored. Recording all points in $k$-space on a uniform grid allows estimation of $\rho$ by inverting a single Fourier transform.

While simple computationally, failing to accommodate for off-resonance effects leads to blurring and image distortion in many types of MR scans \cite{Lai}. Caused by static field inhomogeneities, these off-resonance effects manifest themselves in scans with long readout times as well as scans near air/tissue interfaces where magnetic susceptibility variations are exacerbated. This is of specific importance when imaging for surgical planning near the nasal sinuses, ears, and cerebral cortex \cite{Moerland1995} \cite{Neufeld2005}.

Conversely, directing forming $A$ is infeasible for practical image sizes. To produce a $512 \times 512$ image the (dense) $A$ matrix takes on dimensions $512^2 \times 512^2$. Storing $A$ in 32-bit floating point precision requires 256GB of space, prohibitively large on most current machines. As a result, methods of operator compression based on combinations of Fourier operators have been developed to efficient apply $A$ \cite{Fessler2010}.

Once $A$ may be applied efficiently, iterative methods can be used to reliably reconstruct $\rho$ provided enough $k$-space data has been recorded \cite{Sutton2003}. However, the recently developed field of compressed sensing has shown that filling all of $k$-space increases scan time beyond what is needed for exact image reconstruction when off-resonance effects are ignored \cite{Romberg2007}.

In this work we hope to rectify these two issues, demonstrating a method for image recovery with a generalized operator when $k$-space is undersampled. Furthermore, we introduce to the literature a regularization term based on framelets which reduces the computational costs associated with a standard compressed sensing based reconstruction.



\section{Method}

\subsection{Low rank inhomogeneity correction}
The factor $e^{-z(\mathbf{r}_j)t_i}$ in (\ref{aij}) is problematic as $t$ is variable preventing us from absorbing the exponential into $\rho$ and working with a purely Fourier system matrix. Severable low rank approximations allow us to rapidly apply $A$ with little storage overhead. The general form is 

\begin{equation}
e^{-z(\mathbf{r})t} = \sum_{l=1}^L b_{l}(t_i)c_{l}(\mathbf{r}_j)
\end{equation}

Many methods of choosing $b_l$ and $c_l$ exist \cite{Man} \cite{Fessler2005}. The time segmentation approach \cite{Sutton2003} selects a subset of $L$ points, $\{ \hat{t}_l \} \in [0,T]$ and interpolates about an optimal rate map value, $\overline{\omega}$,

\begin{equation}
b_{il} = b_l(t_i)e^{-\overline{\omega}t_i}
\end{equation}

\begin{equation}
c_{lj} = e^{-(z(\mathbf{r}_j) - \overline{\omega}t_i)\hat{t}_l}
\end{equation}

Application of $A$ proceeds as
\begin{equation}
(Ax)_i = \sum_{l=1} b_{il} \sum_{j=1}^n x_j c_{lj} e^{-2\pi i \mathbf{k}(t_i) \cdot \mathbf{r}_j}.
\end{equation}
In the notation of (\ref{dissigeq}),
\begin{equation}\label{aprox}
A = \sum_{l=1}^L B_l \mathcal{G} C_l 
\end{equation}
where $B_l$ and $C_l$ are diagonal matrices ($B_l = diag (b_{il})$, $C_l = diag (c_{lj})$) and $\mathcal{G}$ is a nonuniform discrete Fourier transform of type 2 \cite{Greengard2004}.

While fast to apply, the modified system matrix obeys uncertainty properties differing from the pure Fourier case.
We recall the concept of a {\it restricted isometry constant}, $\delta_s^A$, of a matrix $A$ for an integer $s$ as the smallest number such that
\begin{equation}
(1-\delta_s^A)\|x\|_{l2}^2 \leq \|Ax\|_{l2}^2 \leq (1+\delta_s^A)\|x\|_{l2}^2
\end{equation}
where $s$ is the number of nonzero entries in $x$. When such a constant exists $A$ is said to satisfy the {\it restricted isometry property} guaranteeing us exact image reconstruction from vastly undersampled $k$-space with overwhelmingly high probability \cite{Donoho2004}. In the case of an undersampled Fourier matrix, much work has been devoted to showing that $\delta_s^{\mathcal{G}}$ is small \cite{Candes2006} \cite{Candes.2006}.

In our corrected system matrix we can infer that each term in the sum (\ref{aprox}) has less than ideal isometric properties by examining
\begin{eqnarray}
\|B_l \mathcal{G} C_l x \|_2 & \leq & \|B_l\| \|\mathcal{G} C_l x\|_2 \\ 
& \leq & \max_i(b_{il}) (1+\delta_s^{\mathcal{G}})\max_j(c_{lj})\| x \|_2
\end{eqnarray}
as well as the analogous lower bound. A crude estimate for the restricted isometry constant for an $L=1$ correction to our inhomogeneity is then
\begin{dmath}
\delta_s^{B_l \mathcal{G} C_l} \leq \min( 1-\min_i(b_{il})\min_j(c_{lj}) + \min_i(b_{il})\min_j(c_{lj})\delta_s^{\mathcal{G}}, \; 
1-\max_i(b_{il})\max_j(c_{lj}) + \max_i(b_{il})\max_j(c_{lj})\delta_s^{\mathcal{G}}).
\end{dmath}
Larger values for $L$ increase the corresponding isometry constants and as a result more samples are required for accurate image recovery than in the pure Fourier case. 


\subsection{Frames and Framelets}

Natural images often have sparse representations within some redundant orthogonal basis $X$ for $\R^n$ []. We call such a system a ``tight frame'' if
\begin{equation}
f = \sum_{g \in X} \left< f,g \right> g.
\end{equation}
as $X$ is overcomplete the redundancy often leads to improved image quality in sparse representations \cite{Daubechies2001}.

Recall that a wavelet system is the collection of dilations and shifts of a finite set $\Psi$
\begin{equation}
X(\Psi) = \left\{ 2^{k/2}\psi(2^kx-j) \: : \: \psi \in \Psi,k,j\in \Z \right\}.
\end{equation}
When $X(\Psi)$ is a tight frame the elements, $\psi$, are called ``framelets''. Stacking the basis elements as row vectors of a matrix $D$ we have that $D^tD = I$. However $DD^t = I$ may not hold. The decomposition $Df$ may be computed efficiently via the discrete wavelet transform. We adopt the piecewise linear B-spline framelet in this work \cite{Cai2010a}. More details on framelets may be found in \cite{Caia} \cite{Daubechies2001} and \cite{Cai2010a}.


\subsection{Sparse recovery via $l1$ minimization}

We reduce scan time by appealing to the theory of compressed sensing, downsampling our Fourier matrix uniformly at random and accounting for the missing data with sparsity promoting regularization. The canonical image reconstruction problem is
\begin{equation}\label{csequ}
\underset{\rho}{\operatorname{argmax}} \: J(\rho)  \text{ subject to } A \rho = s
\end{equation}
where $J$ represents some form of $l1$ regularization.

While MR images are not sparse in the image domain, they are sparse in an appropriately chosen transform domain \cite{Lustig2007}. Representation of images in bases of wavelets \cite{Duarte}, and by extension, framelets \cite{Daubechies2001} yield sparse collections of coefficients. Similarly, reconstruction methods based on the total variation norm \cite{Rudin1992} as well as its nonlocal counterpart \cite{Liang2009} have been shown to accurately reproduce detailed images from sparse frequency data.

In this work we advocate a hybrid of total variation and framelet regularization,
\begin{equation}
J(\rho) = | \nabla \rho| + | D \rho |
\end{equation}
where the first term is a total variation norm and $D$ is the discrete framelet decomposition \cite{Cai2008}. The total variation term enhances edges in our reconstructed image while the inclusion of the framelet term allows us to reconstruct smooth images. It has been found previously that hybrid regularizations often lead to superior image quality \cite{Goldstein2009a}.

Not only does this combined regularization improve image quality, it markedly decreases computational cost when compared with a solitary total variation regularizer as the subproblems we solve in our minimization become better conditioned.


\subsection{Split Bregman iterations for image recovery}

Our image is the solution to the constrained optmiziation
\begin{align}\label{jcseq}
\begin{array}{c}
\underset{\rho}{\operatorname{argmax}} \:  | \nabla \rho| + | D \rho | \\
\text{ subject to } A \rho = s
\end{array}
\end{align}
which we solve with the Split Bregman method of Goldstein et al. \cite{Goldstein2009a}. We begin by converting the constrained optimization (\ref{jcseq}) into a sequence of unconstrained problems via Bregman iteration \cite{Osher} 
\begin{align}\label{gbregman}
\left\{
    \begin{array}{rcl}
\rho^{k+1} &=& \min_\rho | \nabla \rho| + | D \rho | + \frac{\mu}{2} \| A\rho - s^k \|^2 \\
s^{k+1} &=& s^k + s - A\rho^{k+1}
\end{array}
\right.
\end{align}
where the parameter $\mu$ affects the convergence rate and is chosen by the user. Alternately updating $\rho^k$ and $s^k$ produces a sequence $\rho^k \rightarrow \rho$, the solution to the constrained problem.

Updating $s^k$ is straightforward. Minimization of the unconstrained step in (\ref{gbregman}) is done by introducing auxiliary variables, $d_x = \nabla_x \rho$, $d_y = \nabla_y \rho$, and $w = D\rho$, allowing us to rewrite our $\rho^k$ update in the equivalent split form

\begin{align}\label{spilt}
    \rho^{k+1} & \left.\begin{array}{l} = \min_\rho \|(d_x,d_y)\|_2 + | w | + \frac{\mu}{2} \| A\rho - s^k \|^2 \\ \end{array}\right. \\
        & 
\text{subject to} \left\{ \begin{array}{l}
       d_x = \nabla_x \rho  \\
       d_y = \nabla_y \rho  \\
       w = D \rho
    \end{array}\right. 
\end{align}

This constrained optimization is then converted to a sequence of unconstrained problems via a second Bregman iteration. This leads us to the following algorithm

\begin{algorithm}[H]
\caption{Split Bregman iteration for constrained optimization}
\label{sbalgo}

\KwInit{ $\rho^0 = A^t s$, and $d_x^0=d_y^0=w^0=b_x^0=b_y^0=b_w^0=0$}

\While{ $\| A\rho^k - s \|_2^2 > tol$ }{
\For{$i=1 \text{ to } n_{inner}$}{

$\rho^{k+1} = \min_\rho \left\{ \begin{array}{l} \frac{\mu}{2}\|A\rho - s^k\|^2 + \frac{\lambda}{2}\|d_x^k - \nabla_x\rho - b_x^k\|^2 \\ + \frac{\lambda}{2}\|d_x^k - \nabla_x\rho - b_x^k\|^2 + \frac{\gamma}{2}\|w^k - D\rho - b_w^k\|^2 \end{array} \right\}$\\

$d_x^{k+1} = shrink(\nabla_x \rho^{k+1} + b_x^k, 1/\lambda)$\\
$d_y^{k+1} = shrink(\nabla_y \rho^{k+1} + b_y^k, 1/\lambda)$\\
$w^{k+1} = shrink(D \rho^{k+1} + b_w^k, 1/\gamma)$\\
$b_x^{k+1} = b_x^{k} + (\nabla_x \rho^{k+1} - d_x^{k+1})$\\
$b_y^{k+1} = b_y^{k} + (\nabla_y \rho^{k+1} - d_y^{k+1})$\\
$b_w^{k+1} = b_w^{k} + (D \rho^{k+1} - w^{k+1})$\\

}
$s^{k+1} = s^k + s - A\rho^{k+1}$\\
}
\end{algorithm}


Here, the function $shrink$ comes from the wavelet literature[] and is defined as
\begin{equation}
shrink(x,a) = \frac{x}{|x|} \max (|x|-a,0).
\end{equation}
The constants $\lambda$ and $\gamma$ are chosen by the user and affect the convergence rate.

Computationally, the $\rho^{k+1}$ update is the most expensive part of our algorithm by a wide margin. The speed of our image reconstruction is thus largely determined by how fast we can solve this minimization. In the purely Fourier case an analytic solution exists leading to a notably fast algorithm \cite{Goldstein2009a}. We have no such formula for the generalized $A$ we work with and instead rely on iterations of the conjugate gradient method to update $\rho^{k+1}$.

By differentiating with respect to $\rho$ and setting the result to zero we find our $\rho^{k+1}$ update as the solution to
\begin{equation}
(\mu A^tA + \lambda \nabla_x^t \nabla_x + \lambda + \gamma D^tD )\rho^{k+1} = rhs^k
\end{equation}
where
\begin{equation}
rhs^k = \mu A^t s^k + \lambda \nabla_x ^t (d_x^k - b_x) + \lambda \nabla_y^t + \gamma D^t(w - b_w).
\end{equation}
Making use of the identities $\nabla^t \nabla = - \triangle$ and $D^tD = I$, gives the system
\begin{equation}\label{hybridsys}
(\mu A^tA - \lambda \triangle + \gamma I )\rho^{k+1} = rhs^k
\end{equation}
which we solve with conjugate gradient iterations.

A major advantage of our hybrid regularizer is now apparent. Consider the system resulting from a regularization based only on total variation (ie $\gamma =0$),
\begin{equation}\label{tvonly}
(\mu A^tA - \lambda \triangle)\rho^{k+1} = rhs^k.
\end{equation}
Denoting maximal and minimal eigenvalues of the matrix in (\ref{tvonly}) by $\lambda_{\max}$ and $\lambda_{\min}$ we can write the condition number of our system as
\begin{equation}
\kappa = \frac{\lambda_{\max}}{\lambda_{\min}}.
\end{equation}

The addition of the framelet regularizer introduces the $\gamma I$ term in (\ref{hybridsys}) reducing the condition number to
\begin{equation}
\kappa_\gamma  =  \frac{\lambda_{\max} + \gamma}{\lambda_{\min} + \gamma} < \kappa
\end{equation}
notably speeding our updates.

\section{Numerical Results}

Our hybrid reconstruction method is compared against three alternate approaches: uncorrected nonuniform-FFT, field corrected total vairation regularization, and field corrected framelet regularization.

The nonuniform Fourier transforms are computed with the min-max method of \cite{Fessler2003} using $6$ points of interpolation. 

\begin{center}
\begin{figure}[h!]
\includegraphics[width=2in,height=2in]{exact_brain_128.eps}
\includegraphics[width=2in,height=2in]{phantom.eps}
\includegraphics[width=2in,height=2in]{fmap.eps}
\caption{Exact $128 \times 128$ brain image, $64 \times 64$ phantom, and field map used in reconstruction experiments. Field map taken from \cite{Fessler2003}.}
\end{figure}
\end{center}


\begin{center}
\begin{figure}[h!]
\includegraphics[width=2in,height=2in]{phantom_actuals_vs_matrix.eps}
\includegraphics[width=2in,height=2in]{phantom_actuals_vs_bregman.eps}
\linebreak
\includegraphics[width=2in,height=2in]{phantom_residuals_vs_matrix.eps}
\includegraphics[width=2in,height=2in]{phantom_residuals_vs_bregman.eps}
\caption{Reconstruction errors on phantom after $2000$ applications of $A$ at $3$x undersampling. Top row: image domain error, $\log \| \rho - \rho_{exact}\|$, vs number of matrix multiplications and Bregman iteration count. Bottom row: residual error, $\log \| A\rho - s \|$. TV, framelet, and hybrid regularization in red, green, and blue, respectively. Splitting parameter $\mu=1$, $\lambda=0.5$, and $\gamma=25$.}
\end{figure}
\end{center}

\begin{center}
\begin{figure}[h!]
\includegraphics[width=1.2in,height=1.2in]{nufft15_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvonly15_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{frameletonly15_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvframelet15_brain_128.eps}
\linebreak
\includegraphics[width=1.2in,height=1.2in]{nufft_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvonly_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{frameletonly_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvframelet_brain_128.eps}
\linebreak
\includegraphics[width=1.2in,height=1.2in]{nufft35_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvonly35_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{frameletonly35_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvframelet35_brain_128.eps}
\linebreak
\includegraphics[width=1.2in,height=1.2in]{nufft45_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvonly45_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{frameletonly45_brain_128.eps}
\includegraphics[width=1.2in,height=1.2in]{tvframelet45_brain_128.eps}
\caption{Brain image reconstruction. Columns correspond to NUFFT, TV, Framelet, and TV-Framelet regularizations. Rows correspond to downsample factors of 1.5 - 4.5x. All iterative methods were stopped after 2997 applications of $A$.}
\end{figure}
\end{center}

\begin{center}
\begin{tabular}{| p{2cm} || r | r | r | r | r | }
  \hline
  \multicolumn{5}{|c|}{Brain image reconstruction error: $\| \rho_{exact} - \rho \|_2$}\\
  \hline                       
  Downsample factor & NUFFT & TV & Framelet & Hybrid \\ \hline
  1.5x & 1.0954e+04 & 0.3134 & 0.3134 & 0.3621 \\
  2.5x & 1.2032e+04 & 32.0884 & 14.3459 & 9.9995 \\
  3.5x & 1.2369e+04 & 35.0548 & 21.8537 & 18.3263 \\
  4.5x & 1.2536e+04 & 48.4947 & 44.9358 & 43.8486 \\
  \hline  
\end{tabular}
\end{center}


\begin{center}
\begin{figure}[h!]
\includegraphics[width=2in,height=2in]{actual_errors_vs_matrix_mult_brain.eps}
\includegraphics[width=2in,height=2in]{actual_errors_vs_bregman_brain.eps}
\linebreak
\includegraphics[width=2in,height=2in]{residual_errors_vs_matrix_mult_brain.eps}
\includegraphics[width=2in,height=2in]{residual_errors_vs_bregman_brain.eps}
\caption{Brain image reconstruction errors after $3000$ applications of $A$ at $2.5$x undersampling. Top row: image domain error, $\log \| \rho - \rho_{exact}\|$, vs number of matrix multiplications and Bregman iteration count. Bottom row: residual error, $\log \| A\rho - s \|$. TV, framelet, and hybrid regularization in red, green, and blue, respectively. Splitting parameter $\mu=1$, $\lambda=0.5$, and $\gamma=25$.}
\end{figure}
\end{center}






\bibliographystyle{plain}
\bibliography{thebib}

\end{document}
